Quantum molecular dynamics simulations for the nonmetal-to-metal transition in 

fluid helium 
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We have performed quantum molecular dynamics simulations for dense helium to study the 
nonmetal-to-metal transition at high pressures. We present new results for the equation of state 
and the Hugoniot curve in the warm dense matter region. The optical conductivity is calculated via 
the Kubo-Greenwood formula from which the dc conductivity is derived. The nonmetal-to-metal 
transition is identified at about 1 g/cm 3 . We compare with experimental results as well as with 
other theoretical approaches, especially with predictions of chemical models. 

PACS numbers: 05.70.Ce, 52.25.Fi, 52.25.Kn, 52.65.Yy, 64.30.+t 
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Hydrogen and helium are by far the most abundant 
elements in nature. The investigation of their phase dia- 
gram, especially at extreme conditions of pressure and 
temperature, is not only of fundamental interest but 
also an indispensable prerequisite for astrophysics. For 
instance, thermophysical properties of hydrogen-helium 
mixtures determine the structure and evolution of stars, 
White Dwarfs, and Giant Planets [J 0, 0. The detec- 
tion of Jupiter-like planets orbiting neighboring stars Q 
has initiated a renewed interest in planetary physics. All 
planetary models require an input of accurate equation 
of state data in order to solve the hydrostatic equations 
for a rotating compact object. While the limiting cases of 
low and high densities are well understood within chem- 
ical and plasma models, the intermediate region is much 
more complex. There, a nonmetal-to-metal transition 
occurs in both hydrogen and helium at pressures of sev- 
eral megabar and temperatures of few eV which implies a 
strong state-dependence of the interparticle interactions 
and, thus, also of the thermodynamic variables. In this 
paper, we study this warm dense matter region where the 
uncertainties in the equation of state data, both experi- 
mentally and theoretically, are greatest. 

A lot of effort has been done to understand the behav- 
ior of warm dense hydrogen [5, 6]. Although the simplest 
element in the periodic table, the transition of a non- 
conducting molecular liquid to an atomic or plasma-like 
conducting fluid at high pressure is still not fully un- 
derstood @, H, . Helium seems to be a much simpler 
system for the study of the high-pressure behavior since 
no dissociation equilibrium between molecules and atoms 
interferes with the ionization equilibrium, and the first 
(24.6 eV) and second ionization energy (54.4 eV) are well 
separated. Surprisingly, only few experimental and the- 
oretical studies exist for warm dense helium 

Gamma. 



In this paper, we perform the first comprehensive ab 
initio study of the high-pressure behavior of helium. We 
determine the equation of state (EOS) in the warm dense 
matter region by means of quantum molecular dynamics 
(QMD) simulations. The Hugoniot curve is derived and 
compared with experimental points fl3| . other ab initio 
calculation s 11411 . a s well as with results of efficient chemi- 
cal models [lOlllll [l2| . Finally, we calculate the dynamic 
conductivity via the Kubo-Greenwood formula and de- 
rive the static conductivity and compare with shock wave 
experimental data • We locate the nonmetal-to-metal 
transition in the high-pressure phase diagram and discuss 
the related plasma phase transition (PPT). 

QMD simulations are a powerful tool to calculate 
the structural, thermodynamic, and optical properties 
of warm dense matter in a combined first-principles ap- 
proach. Details of this method have been described else- 
where [13, mi m, m . 



We have performed ab initio molecular dynamics simu- 
lations employing a finite temperature Fermi occupation 
of the electronic states using Mermins finite temperature 
density functional theory (FT-DFT) (H. The implemen- 
tation of the QMD method comes from VASP (Vienna 
Ab Initio Simulation Package) , a plane wave density func- 
tional code developed at the University of Vienna (25[. 
All electrons are treated quantum mechanically on the 
level of DFT. The electronic wave function is relaxed at 
each QMD step, which assumes decoupled electronic and 
ionic time scales. We have chosen a simulation box with 
32 to 64 atoms and periodic boundary conditions. The 
electron wave functions are modeled using the projec- 
tor augmented wave (PAW) potentials [26J supplied with 
VASP [25j. These PAW potentials enable more accurate 
results for conductivity calculations compared with other 
pseudopotentials. The exchange correlation functionals 
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are calculated within generalized gradient approximation 
(GGA). Our most accurate calculations were done using 
the GGA parameterization of PBE 27J . The convergence 
of the thermodynamic quantities in QMD simulations is 
of significant importance [28j]. We have chosen a plane 
wave cutoff E cut at 700 eV where the pressure is con- 
verged to within 2%. 

Performing QMD simulations to calculate the EOS of 
He only the T point was used for the representation of 
the Brillouin zone. Calculations for Al have shown that 
it is not recommended to calculate higher-order k-point 
sets [It}- Furthermore, the mean value point (1/4, 1/4, 
1/4) was used for conductivity calculations. 

Our simulations were performed for a canonical ensem- 
ble where the temperature, the volume of the simulation 
box, and the particle number in the box are conserved 
quantities during a QMD run. To keep the tempera- 
ture on a predefined level, the ion temperature is reg- 
ulated by a Nose-Hoover thermostat and the electronic 
temperature is fixed by Fermi weighting the occupation 
of bands (25|. After about hundred time steps the sys- 
tem is equilibrated and the subsequent 400 to 1000 time 
steps are taken to calculate pressures, energies and other 
quantities as averages over this simulation period. 

First we present results for the thermal and caloric 
EOS of warm dense helium in Figs. [T] and [2j The 
isotherms of the pressure and the internal energy be- 
have very systematically with temperature and density 
and show no indications of an instability such as the 
PPT at lower temperatures, contrary to results derived 
within the chemical picture [lCj, [ll|, [l2| . For instance, the 
EOS of Winisdoerffer and Chabrier (WC) [l2j is based 
on a free energy minimization schema for a mixture of 
helium atoms, single and double charged ions, and free 
electrons. Correlations are taken into account based on 
effective two-particle potentials. It agrees well with our 
QMD results for the pressure up to about 1 g/cm 3 and 
for ultra-high densities above about 50 g/cm 3 . How- 
ever, this chemical model shows a systematic trend to 
lower pressures in the intermediate, strongly coupled re- 
gion where the QMD results already approach an almost 
temperature-independent behavior as characteristic of a 
degenerate electron gas. These results underline the sig- 
nificance of ab initio calculations for warm dense matter 
states and will have a strong impact on calculations of 
planetary interiors [29J. Furthermore, we can identify the 
region where efficient chemical models are applicable in 
favor of time-consuming ab initio calculations. 

Based on this EOS data, we have determined the Hugo- 
niot curve via the condition 



(E-E ) = ±(P + P )(V -V), 



(1) 



which relates all final states of a shock wave experi- 
ment (E,P, V) with the initial state (E , P , Vq). In 
our calculations we have used the values Eq = 20 J/g, 
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FIG. 1: Pressure isotherms in comparison with the WC free 
energy model [l2| . 
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FIG. 2: Isotherms of the internal energy. 

T = 4 K, P = 1 bar-C P, and V = 32.4 cm 3 /mol 
(q q = 0.1235 g/cm 3 ) for the first shock and E x = 57 kJ/g, 
Ti = 15.000 K, Pi = 190 kbar, and V x = 9.76 cm 3 /mol 
(qi = 0.41 g/cm 3 ) for the secound shock. We com- 
pare our results with double-shock experiments of Ncl- 
lis et al. [l]|, with recent DFT-MD calculations of Mil- 
itzer 14 1, and with two chemical models 12, 3(3] in Fig. [3] 
A very good agreement is found with the double-shock 
experiments as well as with the other theoretical Hugo- 
niot curves up to about 1 g/cm 3 which results from the 
accordance of the EOS data up to this density as men- 
tioned above. A central problem in this context is the 
value and location of the maximum compression ratio 
?7max- The chemical model FVT^ [3(3] is based on fluid 
variational theory and considers an ideal ionization equi- 
librium in addition. This model predicts a value of 5.5 
at 375 GPa and 50.000 K if the second shock of the ex- 
periment is taken as initial state. Militzer [3] found a 
maximum compression of 5.24 at 360 GPa and 150.000 K 
for the principal Hugoniot starting at (Eq,Pq,Vq) by 
using an EOS composed of zero-Kelvin DFT-MD re- 
sults accounting for excited states for lower tempera- 
tures as shown in Fig. [3] and Path Integral Monte Carlo 
(PIMC) data in the high-temperature limit. Note that 
only a finite-temperature DFT calculation yields the self- 
consistent thermal ground state of the system, which is 
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not equivalent to applying a thermal occupation of the 
empty electronic states ("excited states") obtained in a 
zero-Kelvin electronic structure calculation as performed 
in for details, see [3l| . 

The QMD simulations were performed up to 1.5 g/cm 3 , 
350 GPa, and 60.000 K where the maximum compression 
ratio has not been reached yet. For still higher temper- 
atures the number of bands increases drastically beyond 
the scope of our computer capacity. Besides PIMC simu- 
lations [3] , orbital- free DFT methods may be applicable 
in that high-temperature region 32]. Interestingly, the 
maximum compression ratio for helium (f7 max > 5) is 
greater than that for hydrogen (ry max = 4.25); see e.g. @ 
for a more detailed discussion. 
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FIG. 3: Hugoniot curve for helium: QMD results are com- 
pared with double-shock experiments of Nellis et al. [l3T ] , zero- 
Kelvin DFT-MD results of Militzer [14[| accounting for excited 
states, and two chemical models WC|il| and FVT+ [3(J. 



The dynamic conductivity ak(ui) for one k point is 
derived from the Kubo-Greenwood formula [15j, Il6l 
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where e is the electron charge and m its mass. The sum- 
mations over i and j run over N descrete bands consid- 
ered in the electronic structure calculation for the cubic 
supercell volume fi. The three spatial directions are av- 
eraged by the a sum. F^e^k) describes the occupation of 
the ith band corresponding to to the energy e^k and the 
wavefunction ^k at k. 6 must be broadened because of 
the discrete eigenvalues resulting from the finite simula- 
tion volume 17] . Integration over the Brillouin zone is 
done by sampling special k points fljj ]. 



aM=5> k MFF(k), 



(3) 



where IF(k) is the weighting factor for the respective k 
point. Calculations are usually done at the mean value 
point 0. 
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FIG. 4: QMD results for the static conductivity are com- 
pared with shock-wave experiments of Ternovoi et al. [20| be- 
tween (15 - 25) x 10 3 K and isotherms of the COMPTRA04 
model [341]: temperatures are indicated. 



The dc conductivity follows from Eq. ((2|) in the limit 
lj — > 0. We compare with isentropic compression ex- 
periments of Ternovoi et al. [2p| performed in the range 
(15 — 25) x 10 3 K and predictions of the partially ion- 
ized plasma model COMPTRA04 [H,!!!] in Fig-H The 
experimental points show a very strong increase between 
0.7 and 1.4 g/cm 3 indicating that a nonmetal-to-metal 
transition occurs. Using the Mott criterion for the mini- 
mum metallic conductivity of 20000/£lm also for finite 
temperatures, this transition can be located at about 
1.3 g/cm 3 . The QMD results reproduce the strong in- 
crease found experimentally very well except for the low- 
est density of 0.72 g/cm 3 where the experimental value 
is substantially lower than the QMD result. 

The COMPTRA04 model [H,li| calculates the ioniza- 
tion degree and, simultaneously, the electrical conductiv- 
ity accounting for all scattering processes of free electrons 
in a partially ionized plasma. This approach is able to 
describe the general trends of the electrical conductivity 
with the density and temperature as found experimen- 
tally, see [HI . The isotherms displayed in Fig. [4] cover 
almost the range of the experimental points and agree 
also with the QMD data so that the nonmetal-to-metal 
transiton is described qualitatively very well. 

The strong influence of the temperature on the dc con- 
ductivity in this transition region can be seen by compar- 
ing the QMD results for two temperatures at the same 
density point; see Fig. [4] In order to exclude systematic 
errors from the experimental temperature determination, 
we have performed EOS calculations for the experimen- 
tal points and compare typical values in Table I. A very 
good agreement is obtained so that the discrepancy in the 
conductivity data for the lowest density stems probably 
from the band gap problem of DFT. In order to solve this 
problem, DFT calculations beyond the GGA have to be 
performed by using, e.g., exact exchange formalsims [36| 
or quasi-particle calculations 37J. This is an important 
issue of future work devoted to the nonmetal-to-metal 



4 



transition. 

TABLE I: Comparison of EOS data inferred from the exper- 
iment [2(| and QMD data. 
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The origin of this nonmetal- to- metal transition can be 
elucidated by inspecting the variation of the density of 
states (DOS) along the path of the shock-wave experi- 
ments, see Fig. [5] For the lowest density of 0.72 g/cm 3 , 
a gap still exists in the DOS so that thermally activated 
electron transport occurs as in semiconductors. With in- 
creasing density, the gap region is slightly reduced. The 
main effect is, however, that electronic states fill up the 
region of the Fermi energy with increasing temperature so 
that a higher, metal-like conductivity follows, see also 
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FIG. 5: DOS as function of the energy for three typical situ- 
ations. The Fermi energy is located at zero. 

In summary, we have determined the thermophysical 
properties of dense helium within an ab initio approach. 
The results show clearly the strong influence of quan- 
tum effects and correlations in the warm dense matter re- 
gion. The nonmetal-to-metal transition occurs at about 
1 g/cm 3 , in good agreement with shock wave experimen- 
tal data. These new results will have a strong influence 
on models for planetary interiors. 
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